
##marginal effects plots
## Rueda's data & model (but with default vcovBK confidence intervals added)

# Instantiate
	library(foreign)
	rueda = read.csv("Duplication dataset.csv")
	attach(rueda)

# rename min wage variable
	minwag = rueda$min

library("lmtest")
library("plm")

## Rueda's models

## policies as DV
dvemp4 = plm(govem~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")
dvwelf4 = plm(generos~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")
dvmin4 = plm(minwag~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")

dvemp5 = coeftest(dvemp4, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
dvwelf5 = coeftest(dvwelf4, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
dvmin5 = coeftest(dvmin4, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))

## policies as IV
ivemp4 = plm(ia5010~govem+generos+minwag+hkcorp+govem:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")
ivwelf4 = plm(ia5010~govem+generos+minwag+hkcorp+generos:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")
ivmin4 = plm(ia5010~govem+generos+minwag+hkcorp+minwag:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")

ivemp5 = coeftest(ivemp4, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
ivwelf5 = coeftest(ivwelf4, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
ivmin5 = coeftest(ivmin4, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))


## change PCSEs to Beck-Katz default (diag=F, cluster=c("group","time"))

dvemp2 <- plm(govem~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")
dvwelf2 = plm(generos~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")
dvmin2 = plm(minwag~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")

dvemp3 = coeftest(dvemp2, vcov=function(x) vcovBK(x, diagonal=F,cluster=c("group","time")))
dvwelf3 = coeftest(dvwelf2, vcov=function(x) vcovBK(x, diagonal=F,cluster=c("group","time")))
dvmin3 = coeftest(dvmin2, vcov=function(x) vcovBK(x, diagonal=F,cluster=c("group","time")))


ivemp2 = plm(ia5010~govem+generos+minwag+hkcorp+govem:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")
ivwelf2 = plm(ia5010~govem+generos+minwag+hkcorp+generos:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")
ivmin2 = plm(ia5010~govem+generos+minwag+hkcorp+minwag:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")

ivemp3 = coeftest(ivemp2, vcov=function(x) vcovBK(x, diagonal=F,cluster=c("group","time")))
ivwelf3 = coeftest(ivwelf2, vcov=function(x) vcovBK(x, diagonal=F,cluster=c("group","time")))
ivmin3 = coeftest(ivmin2, vcov=function(x) vcovBK(x, diagonal=F,cluster=c("group","time")))


## marginal effects plots with Rueda & BK confidence intervals

par(mfrow=c(2,3))

## partisan --> policy

## govpart --> govem
plot(na.omit(hkcorp),dvemp2$coeff["govpart"] + na.omit(hkcorp)*dvemp2$coeff["govpart:hkcorp"],"l",lwd=3,ylab="Marginal effect of left government",xlab="Level of corporatism",main="Marginal effect of left govt on govt employment")
## BK conf int
lines(lowess(na.omit(hkcorp), dvemp2$coeff["govpart"] + na.omit(hkcorp)*dvemp2$coeff["govpart:hkcorp"] + 1.96*sqrt( dvemp3["govpart",2]^2 +(na.omit(hkcorp)^2)*dvemp3["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvemp2, diagonal=F,cluster=c("group","time"))[1,"govpart:hkcorp"]) )),col="red",lty=2)
lines(lowess(na.omit(hkcorp), dvemp2$coeff["govpart"] + na.omit(hkcorp)*dvemp2$coeff["govpart:hkcorp"] - 1.96*sqrt( dvemp3["govpart",2]^2 +(na.omit(hkcorp)^2)*dvemp3["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvemp2, diagonal=F,cluster=c("group","time"))[1,"govpart:hkcorp"]) )),col="red",lty=2)
## Rueda conf int
lines(lowess(na.omit(hkcorp), dvemp4$coeff["govpart"] + na.omit(hkcorp)*dvemp4$coeff["govpart:hkcorp"] + 1.96*sqrt( dvemp5["govpart",2]^2 +(na.omit(hkcorp)^2)*dvemp5["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvemp4, diagonal=T,cluster="time")[1,"govpart:hkcorp"]) )),col="red")
lines(lowess(na.omit(hkcorp), dvemp4$coeff["govpart"] + na.omit(hkcorp)*dvemp4$coeff["govpart:hkcorp"] - 1.96*sqrt( dvemp5["govpart",2]^2 +(na.omit(hkcorp)^2)*dvemp5["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvemp4, diagonal=T,cluster="time")[1,"govpart:hkcorp"]) )),col="red")
## show threshold for stat signif of marg effect
abline(h=0,col="gray")
## high & low corp lines
abline(v=.15,col="blue",lwd=2)
abline(v=.9,col="blue",lwd=2)
 
## govpart --> generos
plot(na.omit(hkcorp),dvwelf2$coeff["govpart"] + na.omit(hkcorp)*dvwelf2$coeff["govpart:hkcorp"],"l",lwd=3,ylab="Marginal effect of left government",xlab="Level of corporatism",main="Marginal effect of left govt on welare generosity", ylim=c(-.1,.1))
## BK
lines(lowess(na.omit(hkcorp), dvwelf2$coeff["govpart"] + na.omit(hkcorp)*dvwelf2$coeff["govpart:hkcorp"] + 1.96*sqrt( dvwelf3["govpart",2]^2 +(na.omit(hkcorp)^2)*dvwelf3["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvwelf2, diagonal=F,cluster=c("group","time"))[1,"govpart:hkcorp"]) )),col="red",lty=2)
lines(lowess(na.omit(hkcorp), dvwelf2$coeff["govpart"] + na.omit(hkcorp)*dvwelf2$coeff["govpart:hkcorp"] - 1.96*sqrt( dvwelf3["govpart",2]^2 +(na.omit(hkcorp)^2)*dvwelf3["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvwelf2, diagonal=F,cluster=c("group","time"))[1,"govpart:hkcorp"]) )),col="red",lty=2)
## Rueda
lines(lowess(na.omit(hkcorp), dvwelf4$coeff["govpart"] + na.omit(hkcorp)*dvwelf4$coeff["govpart:hkcorp"] + 1.96*sqrt( dvwelf5["govpart",2]^2 +(na.omit(hkcorp)^2)*dvwelf5["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvwelf4, diagonal=T,cluster="time")[1,"govpart:hkcorp"]) )),col="red")
lines(lowess(na.omit(hkcorp), dvwelf4$coeff["govpart"] + na.omit(hkcorp)*dvwelf4$coeff["govpart:hkcorp"] - 1.96*sqrt( dvwelf5["govpart",2]^2 +(na.omit(hkcorp)^2)*dvwelf5["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvwelf4, diagonal=T,cluster="time")[1,"govpart:hkcorp"]) )),col="red")
## signif
abline(h=0,col="gray")
## high/low
abline(v=.15,col="blue",lwd=2)
abline(v=.9,col="blue",lwd=2)

## govpart--> minwag
plot(na.omit(hkcorp),dvmin2$coeff["govpart"] + na.omit(hkcorp)*dvmin2$coeff["govpart:hkcorp"],"l",lwd=3,ylab="Marginal effect of left government",xlab="Level of corporatism",main="Marginal effect of left govt on minimum wage")
## BK
lines(lowess(na.omit(hkcorp), dvmin2$coeff["govpart"] + na.omit(hkcorp)*dvmin2$coeff["govpart:hkcorp"] + 1.96*sqrt( dvmin3["govpart",2]^2 +(na.omit(hkcorp)^2)*dvmin3["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvmin2, diagonal=F,cluster=c("group","time"))[1,"govpart:hkcorp"]) )),col="red",lty=2)
lines(lowess(na.omit(hkcorp), dvmin2$coeff["govpart"] + na.omit(hkcorp)*dvmin2$coeff["govpart:hkcorp"] - 1.96*sqrt( dvmin3["govpart",2]^2 +(na.omit(hkcorp)^2)*dvmin3["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvmin2, diagonal=F,cluster=c("group","time"))[1,"govpart:hkcorp"]) )),col="red",lty=2)
## Rueda
lines(lowess(na.omit(hkcorp), dvmin4$coeff["govpart"] + na.omit(hkcorp)*dvmin4$coeff["govpart:hkcorp"] + 1.96*sqrt( dvmin5["govpart",2]^2 +(na.omit(hkcorp)^2)*dvmin5["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvmin4, diagonal=T,cluster="time")[1,"govpart:hkcorp"]) )),col="red")
lines(lowess(na.omit(hkcorp), dvmin4$coeff["govpart"] + na.omit(hkcorp)*dvmin4$coeff["govpart:hkcorp"] - 1.96*sqrt( dvmin5["govpart",2]^2 +(na.omit(hkcorp)^2)*dvmin5["govpart:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(dvmin4, diagonal=T,cluster="time")[1,"govpart:hkcorp"]) )),col="red")
## signif
abline(h=0,col="gray")
## high/low
abline(v=.15,col="blue",lwd=2)
abline(v=.9,col="blue",lwd=2)

## policy --> inequality ##

## govem --> inequality
plot(na.omit(hkcorp),ivemp2$coeff["govem"] + na.omit(hkcorp)*ivemp2$coeff["govem:hkcorp"],"l",lwd=3,ylab="Marginal effect of govt employment",xlab="Level of corporatism",main="Marginal effect of govt employment on inequality",ylim=c(-.06,.04))
## BK
lines(lowess(na.omit(hkcorp), ivemp2$coeff["govem"] + na.omit(hkcorp)*ivemp2$coeff["govem:hkcorp"] + 1.96*sqrt( ivemp3["govem",2]^2 +(na.omit(hkcorp)^2)*ivemp3["govem:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivemp2, diagonal=F,cluster=c("group","time"))[1,"govem:hkcorp"]) )),col="red",lty=2)
lines(lowess(na.omit(hkcorp), ivemp2$coeff["govem"] + na.omit(hkcorp)*ivemp2$coeff["govem:hkcorp"] - 1.96*sqrt( ivemp3["govem",2]^2 +(na.omit(hkcorp)^2)*ivemp3["govem:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivemp2, diagonal=F,cluster=c("group","time"))[1,"govem:hkcorp"]) )),col="red",lty=2)
## Rueda
lines(lowess(na.omit(hkcorp), ivemp4$coeff["govem"] + na.omit(hkcorp)*ivemp4$coeff["govem:hkcorp"] + 1.96*sqrt( ivemp5["govem",2]^2 +(na.omit(hkcorp)^2)*ivemp5["govem:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivemp4, diagonal=T,cluster="time")[1,"govem:hkcorp"]) )),col="red")
lines(lowess(na.omit(hkcorp), ivemp4$coeff["govem"] + na.omit(hkcorp)*ivemp4$coeff["govem:hkcorp"] - 1.96*sqrt( ivemp5["govem",2]^2 +(na.omit(hkcorp)^2)*ivemp5["govem:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivemp4, diagonal=T,cluster="time")[1,"govem:hkcorp"]) )),col="red")
## signif
abline(h=0,col="gray")
## high/low
abline(v=.15,col="blue",lwd=2)
abline(v=.9,col="blue",lwd=2)

## generos --> inequality
plot(na.omit(hkcorp),ivwelf2$coeff["generos"] + na.omit(hkcorp)*ivwelf2$coeff["generos:hkcorp"],"l",lwd=3,ylab="Marginal effect of welfare generosity",xlab="Level of corporatism",main="Marginal effect of welfare generosity on inequality",ylim=c(-.007,.007))
## BK
lines(lowess(na.omit(hkcorp), ivwelf2$coeff["generos"] + na.omit(hkcorp)*ivwelf2$coeff["generos:hkcorp"] + 1.96*sqrt( ivwelf3["generos",2]^2 +(na.omit(hkcorp)^2)*ivwelf3["generos:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivwelf2, diagonal=F,cluster=c("group","time"))[1,"generos:hkcorp"]) )),col="red",lty=2)
lines(lowess(na.omit(hkcorp), ivwelf2$coeff["generos"] + na.omit(hkcorp)*ivwelf2$coeff["generos:hkcorp"] - 1.96*sqrt( ivwelf3["generos",2]^2 +(na.omit(hkcorp)^2)*ivwelf3["generos:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivwelf2, diagonal=F,cluster=c("group","time"))[1,"generos:hkcorp"]) )),col="red",lty=2)
## Rueda
lines(lowess(na.omit(hkcorp), ivwelf4$coeff["generos"] + na.omit(hkcorp)*ivwelf4$coeff["generos:hkcorp"] + 1.96*sqrt( ivwelf5["generos",2]^2 +(na.omit(hkcorp)^2)*ivwelf5["generos:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivwelf4, diagonal=T,cluster="time")[1,"generos:hkcorp"]) )),col="red")
lines(lowess(na.omit(hkcorp), ivwelf4$coeff["generos"] + na.omit(hkcorp)*ivwelf4$coeff["generos:hkcorp"] - 1.96*sqrt( ivwelf5["generos",2]^2 +(na.omit(hkcorp)^2)*ivwelf5["generos:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivwelf4, diagonal=T,cluster="time")[1,"generos:hkcorp"]) )),col="red")
## signif
abline(h=0,col="gray")
## high/low
abline(v=.15,col="blue",lwd=2)
abline(v=.9,col="blue",lwd=2)

## minwag --> inequality
plot(na.omit(hkcorp),ivmin2$coeff["minwag"] + na.omit(hkcorp)*ivmin2$coeff["minwag:hkcorp"],"l",lwd=3,ylab="Marginal effect of minimum wage",xlab="Level of corporatism",main="Marginal effect of minimum wage on inequality",ylim=c(-0.6,.4))
## BK
lines(lowess(na.omit(hkcorp), ivmin2$coeff["minwag"] + na.omit(hkcorp)*ivmin2$coeff["minwag:hkcorp"] + 1.96*sqrt( ivmin3["minwag",2]^2 +(na.omit(hkcorp)^2)*ivmin3["minwag:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivmin2, diagonal=F,cluster=c("group","time"))[1,"minwag:hkcorp"]) )),col="red",lty=2)
lines(lowess(na.omit(hkcorp), ivmin2$coeff["minwag"] + na.omit(hkcorp)*ivmin2$coeff["minwag:hkcorp"] - 1.96*sqrt( ivmin3["minwag",2]^2 +(na.omit(hkcorp)^2)*ivmin3["minwag:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivmin2, diagonal=F,cluster=c("group","time"))[1,"minwag:hkcorp"]) )),col="red",lty=2)
## Rueda
lines(lowess(na.omit(hkcorp), ivmin4$coeff["minwag"] + na.omit(hkcorp)*ivmin4$coeff["minwag:hkcorp"] + 1.96*sqrt( ivmin5["minwag",2]^2 +(na.omit(hkcorp)^2)*ivmin5["minwag:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivmin4, diagonal=T,cluster="time")[1,"minwag:hkcorp"]) )),col="red")
lines(lowess(na.omit(hkcorp), ivmin4$coeff["minwag"] + na.omit(hkcorp)*ivmin4$coeff["minwag:hkcorp"] - 1.96*sqrt( ivmin5["minwag",2]^2 +(na.omit(hkcorp)^2)*ivmin5["minwag:hkcorp",2]^2 + 2*na.omit(hkcorp)*( vcovBK(ivmin4, diagonal=T,cluster="time")[1,"minwag:hkcorp"]) )),col="red")
## signif
abline(h=0,col="gray")
## high low
abline(v=.15,col="blue",lwd=2)
abline(v=.9,col="blue",lwd=2)

